Sequential fragmentation: The origin of columnar quasi-hexagonal patterns 
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We present a model that explains the origin and predicts the statistical properties of columnar 
quasi-hexagonal crack patterns, as observed in the columnar jointing of basaltic lava flows. Irregular 
fractures appear at the surface of the material, induced by temperature gradients during cooling. At 
later times fractures penetrate into the material, and tend to form polygonal patterns. We show that 
this ordering can be described as a tendency to minimize an energy functional. Atomistic simulations 
confirm this interpretation. Numerical simulations based on a phenomenological implementation of 
this principle generate patterns that have remarkably good statistical agreement with real ones. 
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I. INTRODUCTION 

Columnar jointing in some kinds of volcanic rocks- 
especially basaltic lava flows-is one spectacular exam- 
ple of geometrical order in nature, where cracks split the 
rock in a set of parallel columns ||^,|^ . Perpendicularly to 
the columns the fractures show a distinctive pattern of 
mostly pentagonal and hexagonal polygons whose sizes 
vary from a few centimeters to about two meters (see 
Fig. It has been realized for more than a century 
that columns result from the contraction of the cooling 
lava after solidification. There is consensus by now that 
fractures start at the surfaces of the igneous body, and 
propagate to the interior as the rock cools down. Frac- 
ture patterns at the surface are rather disordered, but 
the fractures progressively order as they penetrate the 
sample, reaching an almost stable polygonal configura- 
tion after some depth. The fundamental reason of this 
ordering process is unknown at present. Further evidence 
for these facts come from the reproduction of columnar 
cracking in samples of dessicating starches In this 
case the columns have diameters in the range of the mil- 
limeters, and the ordering process is apparent. 

Ryan and Sammis [Q collected evidence suggesting 
that the vertical advance of the fracturing front is not 
continuous, but a discrete process that we call sequential 
fragmentation: at each step, and when some maximum 
tensile stress is exceeded, a layer of material is fractured. 
The existence of this punctuated advance of the fractur- 
ing front can be seen in the lateral faces of the columns, 
which usually have typical marks called striae, or chisel 
marks Further work |^,^ has clarified the form in 
which fractures propagate within each horizontal layer. 
A fracture appears at some point and propagates un- 
der the combination of mode I and mode III fracturing, 
guided by the upper part of the rock, which was frac- 
tured previously. This is a justification for the prismatic 
form of the columns and the striae on their faces. How- 
ever, it assumes that the polygonal pattern is already 
formed in the upper parts of the rock. Aydin and De- 
Graff H tried to give a justification for the evolution 



of superficial fractures, which usually meet in the form 
of a 'T', towards the 'Y' triple junctions typical of well 
developed columns. But, as they point out, the predic- 
tion of the overall polygonal pattern of fractures would 
require a three dimensional mechanical analysis of the 
interaction among many neighbor triple junctions. This 
is a extremely difficult task if pursuit from a microscopic 
point of view. They conclude that probably energetic 
arguments (involving fracture energy and elastic energy) 
may dictate the way in which the final polygonal pattern 
is formed. 

Energetic arguments have been invoked since quite a 
long time to justify the polygonal structure of columnar 
basalts 0, and it is known that the perfect hexagonal 
pattern relieves the maximum amount of elastic energy 
for a given total length of fractures We argue in the 
next section that for the problem of sequential fragmen- 
tation (in which new parts of the rock are sequentially 
fractured under the infiuence of both the previously frac- 
tured parts and the still intact rock underneath) a min- 
imum principle can be invoked to describe the evolution 
in depth of the fracture pattern. In the next section we 
actually show in a numerical simulation (which assumes 
only the sequential nature of the fragmentation process 
and no other ad-hoc assumption) that the fracture pat- 
tern starts being disordered at the surface, and progres- 
sively orders as it penetrates the sample, approaching 
the ideal pattern we expect on an energetic basis. The 
minimum principle is used in section III in phenomeno- 
logical simulations, to evolve superficially disordered pat- 
terns into stable polygonal configurations. The statistical 
properties of the final patterns agree with the available 
experimental data in basalts and also in starches. In sec- 
tion IV we summarize our results, and point out some 
open problems in columnar jointing, mainly associated 
to realistic conditions of cooling. 
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II. ENERGETIC DESCRIPTION OF THE 
ORDERING PROCESS 

We will concentrate on the problem of a semi-infinite 
solid body ('the rock') cooling down through a (horizon- 
tal) free surface ^ . Since this is a situation of inhomo- 
geneous cooling, there will be thermal gradients within 
the rock. Thermal gradient will point vertically at ev- 
ery point, and then temperature will be constant in any 
horizontal plane. Under the stresses generated by the 
thermal gradient, the rock will fracture. 

There are two qualitatively different stages in the frac- 
turing process. One is the appearance of fractures at the 
surface of the rock. Here, the first fracture appears when 
some maximum stress is exceeded at some point of the 
sample. Then it propagates horizontally under the in- 
fiucnce of the inhomogeneities of the rock. When new 
fractures nucleate at the surface, they propagate until 
they meet older ones, usually at right angles, giving rise 
to typical surface fragmentation patterns that have been 
extensively studied, both experimentally [|io|,0 and the- 
oretically [l^ , [r^ . For our purposes, we only mention that 
this stage is governed to a large extent by the random dis- 
order present in the system, since fractures nucleate at 
points where the body can resist the lowest strain. The 
pattern at the surface is usually quite disordered. 

In this paper we study the second stage of the frac- 
turing process of the rock, namely the way in which the 
superficial, disordered pattern of fractures penetrates the 
body and orders. We will assume that the temperature 
distribution within the rock is a given function of coor- 
dinates and time, independent of the actual arrangement 
of fractures, and homogeneous at each horizontal plane. 
The last fact, however, is not enough to assure that the 
fracture front (i.e, the vertical coordinate up to which 
fractures have penetrated, as a function of the horizontal 
coordinate) will be horizontal. In fact, in a standard sit- 
uation of fracture mechanics [Q , it would occur that as 
soon as a fracture penetrates slightly more than the rest, 
stresses accumulate onto that fracture, the result being 
that typically a single fracture advances. For our case 
however, and under realistic cooling conditions, the tem- 
perature gradient decreases ahead of the fractures, so if 
a fracture advances, it rapidly reaches regions where the 
lower temperature gradient precludes the further advance 
of that fracture. This is the reason for the sequential ad- 
vance of the fracture front as schematically illustrated in 
Fig. ^. At each "time step" , a horizontal slab of material 
right below the fracture front is fractured under the in- 
fiuence of the already fractured material above, and the 
still unfractured material below. We will always assume 
that the conditions for sequential fragmentation apply. 

/^From now on wc will treat the rock as a collection 
of particles, elastically joined to their nearest neighbors, 
in the presence of a constant temperature gradient in 
the vertical direction. Changes in temperature are in- 
terpreted as changes in the equilibrium distance between 



particles. Fractures will be modeled by the saturation of 
the elastic energy between neighbor particles as they are 
taken apart a distance larger than some pre-fixed value 
dcr{T) (see Fig. ^). Note that, in this way, for any 
given temperature distribution, and for any arrangement 
of particles, we can define a total energy E for the system 
without ambiguity. 

It is useful to divide the total energy E in two parts, 
E — El + E2. The Ei term, that we call elastic en- 
ergy, comes from those particles being at relative dis- 
tance lower than dcr- This is an elastic energy since 
is quadratic in the relative distance between particles. 
The second part E2 is the contribution from particles 
at relative distance larger than dcr, and then it can be 
associated with broken links (since in this case force van- 
ishes), and identified with the fracture energy. Actually, 
whenever we talk of the existence of a fracture at a given 
position of the sample, we mean that neighbor particles 
are separated by a distance larger than dcr across that 
'fracture'. 

Having defined the degrees of freedom of the system 
and the total energy, we can think of the system as a 
point V in the configuration space of all particle coordi- 
nates. The sequential evolution we have described (Fig. 

corresponds to the sequential mechanical relaxation 
of all particles within the slab between Zi and ^i+i, with 
dz = Zi+i — Zi being the thickness of the slab being frac- 
tured at step i. This sequential process corresponds to 
the movement of V in the energy landscape. Assuming 
that the mechanical relaxation occurs by some kind of 
'viscous' dynamics, the present description becomes com- 
plete and deterministic, and then we can solve in princi- 
ple all the (non-linear) mechanical equations for the prob- 
lem, and obtain in all detail the way in which fractures 
penetrate the sample. The qualitative features of this ad- 
vance, however, can be inferred from general arguments. 
In fact, with the fracture front at a given z position, we 
can calculate the stress field ahead of the fractures, and 
determine the directions along which this stress is max- 
imized. These are the directions that fractures have the 
tendency to follow as they advance. The system releases 
the maximum amount of energy when fractures advance 
along these directions, compared to any other. In other 
words, at each step the configuration point V moves fol- 
lowing a steepest descendent path in the potential energy 
landscape. Note that due to the particular conditions of 
sequential advance, this movement is 'quasistatic', in the 
sense that it does not involve the runaway of fractures 
ahead the fracture front. 

The kind of argument we are using is equivalent to 
those used in surface fragmentation to justify the fact 
that new fractures meet older ones at right angles [||. 
This is a consequence of the tendency of fractures to 
advance perpendicularly to the direction of maximum 
stress, and is equivalent to say that the configuration 
point V moves down in energy following the steepest de- 
scendent path. We are just saying that for sequential 
fragmentation the advance of all fractures is governed 
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by this kind of principle. Then, our minimum principle, 
central to all this work, states that under sequential frag- 
mentation conditions the advance of the fracture front 
occurs with a tendency to reduce as much as possible 
the total energy of the system. Note that during this 
ordering process, the existence of small inhomogeneities 
in the material plays no significant role, as energy will be 
mostly dependent only on the geometrical configuration 
of fractures. 

Our principle then justifies qualitatively the observed 
tendency to produce polygonal arrangements. It is im- 
portant to note, that the system finds the most conve- 
nient pattern by modifying the one at the surface (which 
usually is quite disordered) through small steps as frac- 
tures penetrate the sample. In the next two sub-sections 
we present results that confirm the validity of our inter- 
pretation. 

A. A stress calculation 



to follow as they advance. Starting at the tips of the 
fractures, we see that these directions go away from each 
other, as indicated by the arrows. This indicates that, if 
sequential fragmentation occurs, the close fractures will 
advance with a tendency to separate from each other, 
and eventually to produce a set of evenly spaced frac- 
tures. In fact, only when the evenly spaced configuration 
is reached, the maximum stress direction will coincide 
with the vertical direction, and from here the pattern is 
not modified. This standard calculation coincides quali- 
tatively with what expected from the minimum principle, 
since a set of evenly spaced fractures is the configuration 
that releases the maximum amount of energy (this is the 
equivalent of the honeycomb lattice in three dimensions 
|16| ). Then we see that the conclusions from our min- 
imum principle do not contradict those obtained from 
more standard analysis. The advantage, however, is that 
the minimum principle is much easier to implement in 
cases where a calculation of stresses is not feasible. 



First of all we want to show that standard stresses 
calculations are consistent with the minimum principle. 
We have calculated the stress field surrounding a system 
of unevenly spaced fractures in two dimensions. More 
specifically, we want to calculate the stress field for a set 
of fractures as depicted in the inset of Fig. ^, namely, 
there are pairs of fractures separated by a distance di, 
and the pairs themselves are separated by some other 
distance c?2- 

We start with lattice points joined by springs to form 
a triangular lattice, then modeling an homogeneous and 
isotropic material with Poisson ratio a — 2/(4 + 
\/3) '-^ 0.35. We simulate a piece of size Ix x Iz in the x 
and z directions respectively, taking periodic boundary 
conditions in the x direction, and open boundary condi- 
tions in the z direction. The springs have a rest length 
do that depends on its vertical coordinate z according to 



do{z) = (ioo(l - Py) 



(1) 



In this way we model a constant temperature gradient in 
the z direction (in the simulations we will use /3 = 0.01). 
The periodic boundary conditions in the x direction are 
taken in such a way that the particles at z = are nomi- 
nally at zero strain, whereas all planes on top of that are 
strained with respect to the preferred distance do. The 
two fractures are introduced in the system by eliminating 
all springs that go across the fractures. 

We have solved numerically the problem, by relaxing 
(with a viscous dynamics) the coordinates of the par- 
ticles in order to obtain the equilibrium configuration. 
Then the stress tensor |15 was calculated and diagonal- 
ized at each position. In Fig. |4| we show the results. At 
each point, the tangent to the line shown in that figure 
is the direction perpendicular to the eigenvector corre- 
sponding to the maximum eigenvalue of the stress ten- 
sor, and then it is the direction that fractures will tend 



B. Atomistic simulations of ordering 

The second result presented to validate the minimum 
principle is an atomistic numerical simulation in three 
dimensional systems. We implement sequential fragmen- 
tation in the following way. We use a generalization of 
the procedure extensively used to study surface fragmen- 
tation ||l2| . In that case a layer of material shrinks while 
it is attached to a fixed underlying layer. We take an 
hexagonal plane of particles, with particles attached to 
their neighbors by generalized springs (with an energy- 
displacement relation as that of Fig. ||) of spring constant 
K and initial natural length do. Their positions are the 
dynamical variables. They are attached to an underly- 
ing hexagonal plane of particles (which are kept fixed to 
their original positions during the simulations) by ver- 
tical springs of constant k. The vertical springs do not 
break. Simulation proceeds by reducing the equilibrium 
distance of the horizontal springs of the layer being sim- 
ulated. The first fracture appears when the equilibrium 
distance between two particles becomes grater than the 
corresponding critical distance dcr of the spring that joins 
them. 

For our simulations of sequential fragmentation, the 
only difference is that we consider also the simulated 
plane to be joined to an upper plane of fixed particles 
by spring of constant fc, and that we simulate the frac- 
turing of the system as a sequence of independent two 
dimensional fragmentation processes. In the simulation 
of the successive layers the position of particles in the 
upper plane are taken equal to the final positions of the 
simulation of the previous layer. In the simulation of the 
first layer, we do not have an upper plane. However, to 
avoid introducing disorder into the system (and in or- 
der to break the homogeneity that would occur for an 
absolutely perfect system) we take an upper plane con- 
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sisting of particles located at the hexagonal lattice plus 
some random displacement, independent for each par- 
ticle. We took this displacement to be 0.5 of the lattice 
parameter. We want to mention that other simulations in 
which disorder was included, and the first layer was sim- 
ulated without any upper plane produced qualitatively 
the same results. The equilibrium distance between par- 
ticles do within the layer being simulated is quasistati- 
cally reduced from some initial value doo to pdoo- We use 
p = 0.89. In the energy of the horizontal springs (Fig. ||) 
we use dcr = do + We also take K/k — 100. 

In Fig. H we see the final pattern of fractures for pro- 
gressively deeper layers n, for a system of 1600 particles. 
As we see, the pattern of fractures that appears is highly 
disordered for the first few planes, with many fractures 
ending in the middle of the sample. When we go inside 
the material, there is a clear tendency to order, forming a 
polygonal pattern reminiscent of the experimental obser- 
vations in basalts. Although in Fig. || some influence of 
the hexagonal structure chosen for the underlying lattice 
is observable, we have verified that the same qualitative 
process of ordering is found also for other underlying ge- 
ometries, namely square. We have also looked at the final 
energy the pattern gets after fracturing, and this quan- 
tity is plotted in Fig. ^ as a function of n. As we see, 
this quantity has a tendency to be minimized as succes- 
sive layers are fragmented, which is the right tendency 
predicted by our arguments. Moreover, in Fig. ^ we also 
plot the energy expected for a perfect hexagonal pattern, 
with the size of the hexagons chosen precisely in order 
to minimize the energy. We see that the solution that 
was found by the system was not the perfect one, but 
very close in energy to that. This is a further confirma- 
tion that the tendency to minimize the final energy is in 
fact the driving force for the formation of the polygonal 
pattern. 

III. PHENOMENOLOGICAL CALCULATION 

Having identified the reason why a superficially disor- 
dered pattern shows a tendency to order as it penetrates 
the material does not exhaust the interesting features of 
the problem. Here we will address the observation that 
patterns are usually seen to be polygonal, but not per- 
fectly hexagonal, as it would be preferred by purely ener- 
getic reasons. We will show that this is a consequence of 
the minimization process, since the system is usually not 
able to reach the absolute minimum of the energy poten- 
tial, but gets trapped in a relative minimum. Since the 
problem becomes computationally too costly to be tack- 
led by the methods of the precedent section, we look for 
a phenomenological approach. We will need to calculate 
in some approximate manner the energy of the system as 
fractures advance, in order to search for fracture patterns 
that tend to minimize the energy. 

A realistic calculation is rather complicated and it will 



be presented elsewhere. Here we will restrict to an heuris- 
tic analysis that however is able to show many of the 
known physical properties of fracture patterns. 

Let us suppose that fractures divide the system in sec- 
tors of well defined areas Ai. We are interested in the 
elastic energy Ei of the system after a vertical advance 
dz of the fractures. To lowest order this energy must be 
a function of the Ai^ of the elastic constants, and of the 
precise thermal state of the material. We will use the 
following expression 

i?i = i?o + 7E^«'^^' (2) 

i 

where 7 > and v > 1 are constants, and we have col- 
lected within i^o all possible terms that do not depend 
on Ai. Three main facts have been used in construct- 
ing expression (^). First, the energy is an independent 
sum over different columns of terms that depend on Aj. 
This is the lowest order contribution we expect, in which 
we disregard contributions proportional to the particu- 
lar form of the columns, and interaction terms between 
neighbor columns. Second, the final energy Ei increases 
if Ai increases (i.e., 7 > 0). This is the right tendency, 
since the final elastic energy becomes lower if new frac- 
tures are introduced in the system, and this implies a 
reduction of the typical Ai. Third, the exponent v must 
be greater than one. This condition implies the tendency 
of the system to make the distribution of Ai as uniform 
as possible in order to reduce Ei. With illustrative pur- 
poses, in the rest of this paper we will use v = 2. We 
have repeated the simulations with v in the 1.5 to 2.5 
range with no significant change. The precise properties 
of the material an the thermal state of the system are 
contained in the value of 7. 

Expression (||) for the final elastic energy has to be 
added with the change in the fracture energy E2 during 
the vertical advance. This is simply given in terms of the 
energy needed to create the new fractures as 

5E2 = rjLdz, (3) 

where 77 is the fracture energy per unit area, and L is the 
total length of fractures perpendicular to the propaga- 
tion direction. Collecting the elastic (||) and fracture (||) 
energy terms, we can rephrase the minimum principle in 
the following form. Upon fracture advance, the energy 
functional 

£ = 7^^? + '?^ (4) 

i 

tends to be minimized. The absolute minimum of (Q) is 
attained by a perfect pattern of hexagons of side 

Un - {"^v/dif (5) 

Now, we will use functional (|^) to evolve irregular pat- 
terns (representing superficial fractures) up to point in 
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which they stabiHze, and then compare their statistical 
properties with real ones. Since we are not able to man- 
age a completely general case, we chose a simple possi- 
bility that turns out to produce quite interesting results. 
We generate the pattern at the surface by a process of nu- 
cleation of linear fractures: from randomly chosen points 
within the plane we propagated two opposite, straight 
fractures. The process was repeated many times, with 
new fractures stopping as soon as they reached an older 
fracture. In Fig. 0(a) we show a typical pattern gener- 
ated by this process [ISj . We simulate the modification of 
the pattern with an algorithm that makes small changes 
to the positions of the nodes at which fractures join. Each 
step in the modification of the pattern corresponds to the 
fracture pattern developing into the rock. The new po- 
sition for a node was accepted if the new value of the 
energy, as given by Eq. (^) was lower than the previous 
value. In addition, at each step of the simulation the 
configuration was checked for the existence of very close 
nodes that can allow a change in the topology of the pat- 
tern according to the sketch of Fig. 0(d). Again, the 
changes were accepted only if they reduce the value of 
E. These processes are important since they change the 
number of sides of the polygons, and allow for a progress 
towards more stable patterns. 

An intermediate pattern in the evolution process is 
shown in Fig. 0(b), and the final one (after which all 
proposed changes of the positions of the nodes increase 
the energy) is shown in Fig. |^(c). Since 'time' on our 
simulations corresponds to 'depth' in the rock, the or- 
dering of our patterns represents the progressive order of 
the real lava fractures deeper into the rock |]l9| . The final 
pattern of Fig. |^(c) is not perfectly hexagonal, and thus 
it is only a relative minimum of There is one single 
effective parameter in the simulation, that can be taken 
to be the side of the perfect hexagonal pattern of mini- 
mum energy l-axin- For our simulations this value, as given 
by (^, is indicated in Fig. 0. The qualitative similarity 
of the final pattern with that of the Giant's, Causeway 
shown in Fig. |l|, is apparent. This polygonal pattern is 
now exposed at the surface of the rock, but there is evi- 
dence that this is not the original surface. In Fig. ^ we 
show two quantities that are a measure of the statistical 
similarity between our patterns and the real ones. In Fig. 
^(a) we see the results for the frequency of appearance 
of polygons with a given number of sides (in this case we 
also include the results on cornstarch by Miiller ) , and 
in Fig. ^(b) the corresponding values for the mean area 
of polygons with a given number of sides, both in our 
simulations and in the real patterns. The configurations 
generated by our model are remarkably realistic. We see 
that, both in real cases and in our simulations, the frac- 
tures never reach a perfect hexagonal pattern. Instead, a 
reproducible distribution of polygons, most of them with 
5, 6, and 7 sides is obtained, with a minor contribution of 
polygons with 4 and 8 sides. Also, polygons with higher 
number of sides have larger area as Fig. Wh) shows. 



IV. SUMMARY AND PERSPECTIVES 

In this paper we have given a first approach to a con- 
sistent model for the existence of columnar polygonal 
patterns in lava flows and some dessicating materials. 
We have shown in numerical simulations on a discrete 
model that fractures appear as irregular cracks at the 
free surface of the material and become ordered as they 
penetrate into the interior. We have argued that this 
effect is a consequence of a tendency to minimize an en- 
ergy functional. The process of minimization follows a 
rough landscape, and is always towards a local minimun. 
This process is therefore monotonically decreasing with 
the thermal fluctuations playing no role. Relying on this 
principle, we showed that the statistical properties of ex- 
perimental polygonal patterns can be reproduced. 

There are still some problems that deserve further con- 
sideration and that we plan to discuss in a forthcoming 
publication. They have to do mainly with the realis- 
tic conditions of cooling. As discussed in section II, it 
is precisely the decreasing of the temperature gradient 
ahead of the fractures that makes possible the sequential 
advance of the fracture front, in a coordinated way all 
across the sample. The detailed study of this problem 
provides predictions for the width of chisel marks on the 
columns. 

We also have to determine in a realistic situation the 
value of the constant 7 and u in expression |[ This will 
allow to calculate in particular the typical values for the 
polygons in basalts and starches. Under realistic cooling 
conditions we also have to face the problem that tem- 
perature changes with time, and the effect of this on the 
advance of the fracture front has to be discussed. 
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FIG. 1. Polygonal pattern seen perpendicularly to some of 
the columns of the Giant's Causeway, a Tertiary lava flow in 
Antrim, Northern Ireland, from Ref. (originally from a 
map by O'Reilly H)). 
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FIG. 2. Schematic representation of sequential fragmenta- 
tion in a two-dimensional geometry . Continuous lines rep- 
resent the advancing fractures. At each time step a layer of 
material of thickness dz fractures. New fractures appear be- 
low those already present, but slight modifications in their 
positions are possible, and in fact crucial to the ordering pro- 
cess. 
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distance 

FIG. 3. Schematic energy vs. distance curve for neighbor 
particles in the discrete model of the rock. There is an equilib- 
rium distance between particles that depends on temperature 
do{T). Deviations from this distance costs an clastic energy 
which is quadratic in the displacement. If the distance be- 
comes greater than some critical distance dcr{T) then energy 
saturates. In this way we model fractures, since in this range 
there is no force between particles. 




FIG. 4. Stresses ahead of an unevenly set of fractures as 
depicted in the inset. In the main figure we plot the stress 
field in the neighborhood of a couple of close fractures (see 
text for details). Directions of maximum stress at the tips of 
the fractures are indicated by the arrows. The simulated box 
is marked in gray in the inset. Periodic boundary conditions 
are used along x, and free boundary conditions along z. The 
dotted box in the inset is the region plotted in the main figure. 
We have used di/di =4. 




FIG. 5. The final patterns of fractures for progressively 
deeper layers n in a sequential fragmentation process for lay- 
ers with 1600 particles. For clarity reasons, in the plots the 
system has been duplicated both in horizontal and vertical di- 
rection (periodic boundary conditions are used). In the plots, 
each thick fracture is formed by small lines mostly perpen- 
dicular to the fracture that join the ends of springs that have 
failed after a contraction up to 0.89 of the original distance 
between lattice points. 
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FIG. 6. Energy per particle of the pattern obtained in the 

simulations shown in the preceding figure, as a function of 
the layer index. The ideal minimum of a size-optimized hon- 
eycomb lattice is also shown. Energy is given in units of 
W-'^RdL. 
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FIG. 7. Numerically evolved patterns of fractures. We see 
the original pattern (a), the final (stable) one (c), and one 
intermediate configuration (b). To avoid spurious edge effects, 
only the central region of a simulation performed on a larger 
sample is shown. The numerical algorithm is described in the 
text. Below the final pattern (c), the side of hexagons Imin 
in the expected ideal honeycomb lattice is indicated. In (d) 
we see the kind of processes that allow for a change in the 
number of sides of adjacent polygons. 
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FIG. 8. (a) Histogram for the relative frequency (normal- 
ized to one) of appearance of polygons with different number 
of sides for the Giant's Causeway, columns in cornstarch Q, 
and from our simulations (an average over 10 final configura- 
tions as that of Fig. ^c) is shown), (b) Areas of the polygons 
normalized to the average area for polygons with different 
number of sides (data on cornstarch are not available). 
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